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I. Introduction 

Matter comes in different phases and usually one can 
switch between them by changing the temperature. Close 
to zero temperature thermal fluctuations disappear and 
quantum fluctuations dominate. In this case, by chang¬ 
ing a corresponding control parameter one can induce 
quantum phase transitions (QPTs) between the differ¬ 
ent ground states of quantum systems. QPTs occur in 
many different physical systems, and they attract a lot 
of attention in condensed-matter physics [1]. 

The reason for the recent surge of interest into QPTs 
are new and exotic quantum phases and critical points, 
which cannot be described within Landau’s theory of 
phase transitions, i.e. they cannot be characterized by 
an order parameter. Examples are topologically ordered 
phases [2], quantum spin liquids [3], or deconfined quan¬ 
tum critical points [4, 5]. 

At critical points different parts of the system are 
quantum mechanically strongly correlated and various 
correlation functions show singular behaviour [6]. Sev¬ 
eral years ago, in quantum information theory it was 
suggested that quantum phases and QPTs can be char¬ 
acterized and distinguished in terms of quantum entan¬ 
glement [7, 8]. 

In the present paper we study how different entangle¬ 
ment measures characterize ground state phases within 
simple one dimensional and two dimensional spin mod¬ 
els. We generate these ground states using tensor net¬ 
work techniques. Special emphasis is paid to properties 
of the entanglement monogamy as expressed through the 
Coffman-Kundu-Wootters (CKW) [9] inequality or simi¬ 
lar inequalities. 

Bipartite entanglement is the most studied type of en¬ 
tanglement in quantum information theory [10], however 
even bipartite entanglement measures are under active 
development, especially with respect to the ‘identihca- 
tion power’ of exotic quantum phases. Recently, the con¬ 
cept of ‘entanglement spectrum’ was introduced [11] and 
is now studied intensively [12, 13]. 

Numerous different measures have been proposed to 


quantify entanglement [14]. Those which are studied 
here [6, 15-18] are listed in the Appendix. Studies of 
the entanglement properties of the inany-body systems 
mainly use the entanglement entropy, the one-tangle, the 
concurrence, and the fidelity [6]. There are also inves¬ 
tigations of the multipartite entanglement properties of 
the states (e.g., tripartite entanglement [19] and global 
entanglement [18]). 

It was found in previous studies [6] that entanglement 
measures are able to determine critical properties of the 
systems, in particular the positions of the critical points. 
Recent studies show that it is also possible to extract 
the critical exponents, as was e.g. done using hnite size 
scaling of the Schmidt gap [20] . The techniques to extract 
critical exponents from different entanglement measures 
are still under development. 

In order to simulate the ground states of 1-dimensional 
and 2-dimensional spin models we use a tensor network 
(TN) approach [21, 22]. For a recent review see Ref. [23]. 
The basic idea of TN methods is to represent the wave 
function of a many-body quantum system by a network 
of interconnected tensors. Experience shows that the TN 
class of numerical methods is rather flexible [23]. TNs can 
handle systems in different dimensions, of finite and infi¬ 
nite size, with different boundary conditions and symme¬ 
tries. They are able to model systems of bosons, fermions 
or frustrated spins and can address different types of 
phase transitions. The market of TN methods now pro¬ 
vides several tens of different methods [23], and each of 
them has its own advantages, disadvantages and areas of 
applicability. 

Matrix product states (MPS) are the most famous 
among the TN states [22]. Powerful algorithms such 
as the Density Matrix Renormalization Group (DMRG) 
[24] or Time-Evolving Block Decimation (TEBD) [25] 
can be formulated in terms of MPS. The two-dimensional 
generalization of the matrix product states is called pro¬ 
jected entangled pair states (PEPS) [26]. Details about 
the PEPS and the MPS can be found in Refs. [27-30]. 

There is a large variety of TN methods available in 
order to determine the quantities which characterize a 
quantum state (order parameters, critical exponents, en- 
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tanglement measures). Two questions arise: (1) which 
numerical method is best suited for the simulation of the 
ground state of a particular model? (2) which quantity 
is most efficiently calculated in order to characterize this 
ground state? 

In the present paper we aim to give some input for 
the answer of these questions using the quantum Ising 
model in a transverse field and XXZ models in ID and 
2D geometries. We use numerical methods which are able 
to treat models in the thermodynamic limit: MPS [30] in 
ID and TERG [31], CTMRG [32, 33] for PEPS in 2D. We 
compare our results with results from other studies based 
on other techniques [15-17]. More specifically, we use 
imaginary-time evolution (the one-site update [34, 35], 
TEBD [36] in ID and the ‘simple update’ scheme [28] in 
2D) to find the approximate ground states for the models 
under investigation. Exploiting translational invariance 
of the ground states enables algorithms with reasonable 
requirements for computational resources [31, 37]. 

As for the second question we calculate a set of dif¬ 
ferent entanglement measures such as one-site entangle¬ 
ment entropy, one-tangle, concurrence of formation and 
assistance, bounds on localizable entanglement, local en¬ 
tanglement, negativity and entanglement per bond and 
compare their ‘characterization power’ of the state of the 
system. Moreover, we will compare the entanglement 
properties of the chosen models in ID and 2D geometries. 
In this way we obtain information on the ‘monogamy of 
entanglement’ or entanglement distribution [9]. The en¬ 
tanglement monogamy properties will be studied in some 
detail. 

This paper is organized as follows. In Sec. II, the one- 
site and two-site reduced density matrices (needed for the 
entanglement measures calculation) are determined from 
the MPS and PEPS representations of the ground states. 
Sec. Ill presents numerical results and their interpreta¬ 
tions for quantum Ising and XXZ models. Gonclusions 
are made in Sec. IV. Various entanglement measures are 
briefly listed and discussed in the Appendix. 


II. Translationally invariant tensor network 
methods 


In this section we briefly describe the different ten¬ 
sor network methods we use to obtain the translationally 
invariant ground state wave functions for various spin 
models. We describe the renormalization steps one needs 
to take in order to prevent exponential increase of the 
bond dimensions of the tensors for 2D systems. Fur¬ 
thermore, we show how reduced density matrices are de¬ 
termined using the tensor entanglement renormalization 
group (TERG) or the corner transfer matrix renormaliza¬ 
tion group (GTMRG). We also briefly discuss a renormal¬ 
ization technique for translationally invariant ID systems 
for comparison. All these methods are closely related 
conceptually, but differ in many technical details. 


A. Imaginary-time evolution 


Imaginary-time evolution is one method for finding 
ground-state wave functions [38] . It evolves an arbitrary 
state \4>) (which contains the ground state as a compo¬ 
nent) into the ground state \tfjo) of the Hamiltonian H, 


l^o) = 


lim 

r—¥oo 


exp{—TH)\ijj) 

||exp(-TiJ)|V’)]r 


( 1 ) 


Strictly speaking, this is correct only if the ground state 
is non-degenerate. If the ground state is degenerate one 
obtains an (arbitrary) linear combination of the degen¬ 
erate ground states. Since the imaginary time evolution 
operator is not unitary, normalization must be explicitly 
ensured through the denominator in the above expres¬ 
sion. 

We will first describe imaginary-time evolution of the 
translation invariant PEPS in two dimensions. In both 
dimensions we use periodic boundary conditions. The 
initial random wave function lip) is constructed from a 
product of equal rank-5 tensors A at the lattice sites 
hj, ■ ■ ; 

l'^) = ' • ') | ' (2) 

The tensors contain random entries. The size of each 
physical (spin) index a is two, since we consider spin- 
1/2 systems only. The size of each virtual bond k, I, m, n 
is D. The tensor trace tTr includes summation over all 
spin configurations and over all bond indices. The tensor 
network describing this state is graphically represented 
in Fig. 1. 

We consider spin systems described by Hamiltonians 
with nearest neighbor interactions only. Therefore, the 
Hamiltonian may be decomposed into four terms. 


H = Hk+Hi+Hm + H„. (3) 


These terms correspond to the four bond types shown in 
Fig. 1. Single-site operators in the Hamiltonian may be 
easily incorporated into these terms. The nearest neigh¬ 
bor interactions within each of the four terms commute 
with each other. However, the four terms i7fc, Hm, Hi 
and Hn do not commute. As a consequence, we cannot 
write the evolution operator = exp{—TH) as a prod¬ 
uct of two-site operators. 

However, we may use a first-order Trotter-Suzuki ex¬ 
pansion [39] 

with a ‘small’ imaginary time step Ar in order to write 
Pr approximately as a product of two-site operators 

n . (5) 

time steps 
sites 


where the /i)/ represent the nearest neighbor interactions 
of type A between site i and j {H\ = ). Since 
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FIG. 1. Graphical representation of a translationally invariant 
PEPS. All tensors A are equal. The letters denote the different 
bond types. For a finite lattice the open ends are connected: left 
with right and top with bottom. The octagon encloses the tensors 
M or M defined in Eqs. (6) and (7), respectively. 


imaginary-time evolution is a projective method, Trot¬ 
ter errors do not accumulate during the evolution [40]. 
Higher orders of the Trotter-Suzuki expansion may be 
used in order to achieve a better of convergence. 

In practice, we repeatedly apply imaginary time step 
operators for different ‘directions’ fc, Z, m, or n to a ran¬ 
dom state until convergence is achieved. Convergence is 
judged using suitable criteria, At and the total number 
of time steps N are chosen in order to achieve the de¬ 
sired accuracy. We start from At = 0.1 as the initial 
time step and then progressively reduce it until the state 
does not change any more within specified limits. We 
take this converged state as our approximation for the 
ground state jV'o)- 


B. Update schemes 

After each imaginary time step the size of the tensors 
describing the state increases. Let us consider the evo¬ 
lution of the state over a k bond. Assuming that the 
tensors and correspond to neighbor¬ 

ing sites connected by a /c-bond, the two-site tensor M 
becomes 


M ‘ ^ 

limiriiljmj rij 

Applying the two-site operator to M pro¬ 

duces the ‘evolved’ two-site tensor 


D 

/ ^ ’^klimini'^kljmjrij 


( 6 ) 


II ^ II 

^'^liTniniljmjrij / ^ \Pk)(7i(7j ^^^l^rniniljmjTij 

(Tiaj 


( 7 ) 


of the same size as M. However, reconstructing from 
this tensor the PEPS tensors A using a singular value 
decomposition 


In order to prevent an exponential growth of the tensor 
size during the imaginary-time evolution, the bond size 
must be suitably reduced at each time step. This should 
be done in such a way that the difference between the 
evolved state j')/;®™*-) with increased dimension and its 
approximation with reduced dimension K = 

111 ^evo\.^ _ I ^approx. ^ 11 ^ jg minimal. An imaginary time step 
together with the necessary reduction of the tensor size 
is called an ‘update step’, and the corresponding method 
to reduce the bond size is called ‘update scheme’. 

Different types of update schemes exist. In general, 
in order to implement an update step one needs to take 
into account the whole environment of two evolving ten¬ 
sors [23]. Update schemes that act this way are called 
‘full updates’. They are numerically rather costly. Less 
demanding on the computational resources is the ‘simple 
update’ scheme [28], which takes the environment into 
account approximately; it is based on a generalization of 
a method developed for ID systems called ‘time-evolving 
block decimation’ (TEBD) [25]. There is also a cluster 
update scheme [41], which compromises between the two 
previously discussed schemes. 

In the present work we use the ‘simple update’ scheme 
following Refs. [28, 40]. Here, we would like to comment 
on two important aspects of this algorithm. 

In the ‘simple update’ scheme one introduces bond vec¬ 
tors Afc, Xi, Am, An in addition to the tensors A of the 
standard PEPS. The connection between the tensors in¬ 
troduced in (2) (now we denote these tensors by A^^^) 
and the tensors of new representation A^^^ is given by 

A graphical represen- 
tation of the translationally invariant PEPS with addi¬ 
tional bond vectors is shown in Fig. 2. The small circles 
in the figure correspond to a = The new 

PEPS representation directly corresponds to the canoni¬ 
cal form of MPS introduced by Vidal [25]. 



FIG. 2. Graphical representation of a translationally invari¬ 
ant PEPS with bond vectors. Small circles correspond to \/AT, 
a = k, 1,771,71. The octagon defines the object to be renormalized. 
The extra VAq taken from the external tensors are used to mimic 
renormalization effects due to the environment. 


dxD^ 

k=l 

( 8 ) 


While the introduction of bond vectors appears to be 
a trivial redefinition of the local tensors, it is their role in 
the renormalization or update scheme, which will prove 
to be non-trivial; naively, one would assume that it is M 
defined in Eq. (7) which is renormalized such that the 


increases the size of the fc-bond from D to d x D^. 
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size of the PEPS tensors do not grow. However, it is the 
tensor S defined by 


Si 


minilimiTii 


— \J'J^rrii \J Arii > 


M, 


(9) 


Uminilnmn 


which is renormalized. Here, M is defined as in Eq. (7) 
but M includes all factors yfX necessary due to the re¬ 
definition of the A tensors, 


M, 




— 's/ ^rrii \l^rii ^ 


D 

A, 

k=l 



It is important to note that S contains extra factors of 
•\/A taken from the environment as indicated by the oc¬ 
tagon in the graphical representation of the tensor S in 
Fig. 2. Such an approach is called ‘mean field approx¬ 
imation’ of the environment. However, this statement 
is intuitive and lacks mathematical rigour. The proce¬ 
dure is justified by numerical success. Vidal has proved 
a number of statements which justify this procedure in 
ID for TEBD [25]. 

There is another aspect of the simple update scheme 
that requires comment. Originally the ‘simple update’ 
method was implemented for studying quantum models 
on a honeycomb lattice [28] , where bipartition of the lat¬ 
tice is necessary, i.e. the ground state in PEPS form is 
given by two tensors H, B representing two sublattices. 

One would expect that on a square lattice with transla¬ 
tional invariance the ground state must be translationally 
invariant too, i.e. no bipartition is expected. However, 
in various papers [27, 32, 33] which use ‘simple update’ 
to simulate the ground state on a square lattice, a bi¬ 
partition is introduced even for translationally invariant 
models. It is stated in Ref. [27] that imaginary-time evo¬ 
lution breaks translational invariance of the lattice. The 
motivation of the bipartition is not explained clearly in 
the literature, thus we want to highlight this aspect of 
the ‘simple update’ scheme. 

Our experience shows that, if we start our evolution 
with a random translationally invariant state given by 
random tensor C sitting on each site, the ‘simple up¬ 
date’ leads to a translationally invariant ground state, 
too. The equality of A and B tensors at the end of the 
imaginary-time evolution is very dependent on the con¬ 
vergence criteria and on the bond size D of the PEPS ten¬ 
sors. The resulting A tensors in this case are rotationally 
symmetric with high accuracy (rotation corresponds to 
cyclic permutation of virtual bond indices). Rotational 
symmetry is also ensured by the approximate equality 
of four bond vectors. The accuracy of their equality is 
given by the convergence condition. Note, that for higher 
D it is much harder numerically to obtain approximately 
equal A and B. 

We suggest that the resulting translational invariance 
could be highly dependent on the numerical implementa¬ 
tion of the singular value decomposition procedure used 


during imaginary-time update. In practise, due to a 
gauge freedom the simple update scheme can lead to a 
bipartitioning of the lattice in general, i.e. translational 
invariance would be superficially broken. In fact, trans¬ 
lational invariance is maintained and could be restored 
explicitly using an appropriate transformation. 

The gauge freedom can be easily demonstrated for a 
product of two equal matrices, 

M = AA = {AA){A-^A) = CD. (10) 

The SVD which is used within the simple update scheme 
as indicated in Eq. (8) 

M = UAV'< = {uVA){VAV^) = CD. (11) 

leads to the purely numerical bipartitioning of the tensors 
on the lattice. 

Moreover, we found that the ‘simple update’ can dis¬ 
tinguish ferromagnetic and antiferromagnetic order in 
the state. This order is defined by the sign of the cou¬ 
pling constant in the two-site Hamiltonian that is used 
during the evolution. Thus, the usage of antiferromag¬ 
netic coupling constant will lead to A and B tensors that 
differ with respect to spin-flip transformation. 

Therefore, as an output of the simple update scheme 
for translationally invariant models we obtain the PEPS 
given by one rank-5 tensor and just one unique 

bond vector A = Afe = A; = Am = A„. After 
the completion of the imaginary time evolution, we 
multiply the bond vectors into the tensors, A^^^ = 
^kimn Am , since the bond vectors are not 

needed any more. The resulting tensor network which 
will be used for further tensor contraction algorithms is 
the same as shown in the Fig. 1, however, with new A 
tensors sitting on each site. 

If the ground state is not purely translationally invari¬ 
ant, i.e. has the antiferromagnetic order, PEPS repre¬ 
sentation would then require two tensors A and B to de¬ 
scribe the state. Here, after the completion of the imag¬ 
inary time evolution, we multiply the bond vectors into 
the tensors A and B again. Thus, in the case of lattice 
bipartition the sublattices corresponding to these tensors 
are denoted as A and B, respectively. The resulting ten¬ 
sor network which will be used from now on is shown in 
Fig. 3. 

In the following we will denote tensors obtained from 
imaginary-time evolution and with incorporated bond 
vectors Aa by A and B without bars for simplicity. 


C. Reduced density matrices for 2D systems 

In this subsection we briefly present two different meth¬ 
ods for the calculation of the n-spin reduced density ma¬ 
trices Pn for 2D systems (n < 4). From the reduced 
density matrices we obtain the expectation value of an n- 
spin operator in the standard way: (On) = Tr(0„pn)- 
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FIG. 3. After the imaginary-time evolution in the case of not 
purely translationally invariant state (i.e. describing antiferromag¬ 
netic state) the tensor network assumes a bipartitioned structure 
in terms of tensors A and B. 

The calculation of the density matrices in a tensor net¬ 
work approach involves a tensor trace, the calculation 
of which is exponentially hard in 2D and, therefore, re¬ 
quires renormalization methods (in contrast, for ID sys¬ 
tems the calculation of the reduced density matrices can 
be achieved in polynomial time). The methods we discuss 
here are the tensor-entanglement renormalization group 
(TERG) [31] and the corner transfer matrix renormaliza¬ 
tion group (CTMRG) [32, 33]. 

The essentials of these methods are described e.g. in 
the papers cited above for the calculation of expectation 
values. Here, we present these methods for the calcula¬ 
tion of reduced density matrices. 


The tensors and are located at sites of sublattice 
A and tensors and at sites of sublattice B. For 
simplicity, from now on we will omit the overbars for the 
indices labeling the various T tensors and just keep in 
mind that the virtual indices have dimensions and 
the physical index has the dimension cP. 



FIG. 4. Tensor network and impurity tensors used in TERG. The 
tensor network is in principle infinitely large. The open lines at the 
impurity sites A, B, C, D in the boxed center of the figure indicate 
the physical spin indices of the impurity tensors. The open lines at 
the boundary of the figure are connected to tensors not shown. 


1. Tensor-entanglement renormalization group 


TERG is based on the tensor renormalization group 
(TRG) method introduced by Levin and Nave [37] for 
classical systems. It was modified for quantum systems 
in Ref. [31] using the concept of ‘impurity’ tensors. In 
the present paper, we name ‘impurity’ positions in a ten¬ 
sor network those positions at which spin operators are 
attached or where the physical indices of the network are 
explicitly kept. At all other positions the physical indices 
are summed over. At each site, which is not an impurity 
site, we define the following tensors (see Fig. 4) 

T-- = 4'^* A'^ 

klmn / , k’I'm’n’^klmn’ 


rpb _ _ \ ' R^ 

^ fhnkl / , k'I'm'n'^klmn- 
a 


with the virtual bonds k = k'k, I = I'l, fh = m'm, n = 
n'n. Each index has dimension D^. 

Furthermore, at the impurity sites we define four ‘im¬ 
purity’ tensors T^, , and with physical bond 

a = a'a oi dimension dP as illustrated in Fig. 4, 


(T 

(T 

(r 


)f. _ = (A) 

' kLmn k / 


A\q 
klmn 
B\a 
fnnkl 

C\q 

klmn 
nD\a 


(T * 

k'l'm'n' 


{A) 


klmn 1 




mnkl ’ 


)f. _ = (H) 

'klmn k / 


k'l'm'n' 


{A)l 


(13) 


klmn ’ 


iTn^muM=iBrjn'k>ABrmnkl- 


Now, depending on the sublattice, we perform one of 
the following singular value decompositions (the arrow 
indicates a reshaping of indices). 


Dc 


Ti 


klmn ^ ^ 


1 

mnoc 1 


Da 

^mnkl ^ ^^[nk){lm) ^ ^ ^ ^nkjH^lmjti' 


(14) 


3=1 


The 5* tensors are obtained from the U and Rl ten¬ 
sors of the SVD M = J2U ARl by multiplication with 
•\/A. These decompositions are illustrated in the Fig. 5. 
Analogous SVDs are performed for tensors T^, T®, 



FIG. 5. SVD of the tensors T“ and T^. The indices a and /3 will 
be truncated in order to prevent exponential growth of the tensors 
T defined in Eq. (16) 
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'Y'D 


e.g. 


-Dc 

(rpA\{^'A^A) j.rA , \ ' cAS oAl 

)klmn ' ^^^{t7'^kl){<TA'mn) ^ . '^a^^kla'^aA'mna^ 

a — 1 


(rTiB\(<^B<^B) 

)mnkl ' {aB'nk)(a'^lm) 


Dc 

. qB 2 aBi 

' / ^ ^(7Bnkj3D(j'^lm0' 


3=1 


(15) 


The last step of the TERG procedure is coarse- 
graining, that is the contraction of four S tensors into 
one T tensor 


^ ^ ^nka^mnfU^lm'y^klS' 
klmn 


as illustrated in Fig. 6 . 



FIG. 6. Coarse graining: Contraction of four S tensors. The 
number of T tensors is reduced by a factor of 2. 


Renormalized T^, tensors are determined 

in the same way, e.g. 


_ \ '' qB2 q1 q 4 qA3 

/a^-yS / V 

klmn 


(17) 


using the decompositions (14) and (15). 

From the Eqs. (14) and (15) we realize that the size of 
the virtual bonds a and /3 is D'^ and D'^cP, respectively, 
so that the virtual bonds of the tensors T would grow ex¬ 
ponentially without suitable truncation. In order to pre¬ 
vent the exponential growth we truncate these indices to 
Dc, that is we neglect small singular values in the expan¬ 
sion Eq. (14) and (15); Dc has to be chosen large enough 
to maintain the relevant physical information but small 
enough to stay within acceptable numerical cost. An ac¬ 
ceptable choice for Dc depends on the virtual dimension 
D of the PEPS. In our calculations we use Dc between 
16 and 64. 

After a sufficient number of the TERG transforma¬ 
tions as described above the tensors T^, 
contain all relevant information necessary to calculate 
observables, e.g. the four-site reduced density matrix 




(18) 


which has to be normalized such that Trp = 1. The trace 
includes summations over virtual indices. Two-site and 


one-site reduced density matrices are then easily obtained 
from p by a partial trace. 

TERG transformations are applied until a convergence 
condition is satisfied. After each TERG transforma¬ 
tion we calculate p until we find < e, 

where p^P denotes the reduced density matrix p at the 
n-th recursion step. The matrix norm is implemented as 

ll^ll = ^ij- III practice, we take e between 10 “® 

and 10 ~^°. 

From the two-site and single-site reduced density ma¬ 
trices we calculate the desired physical quantities in sec¬ 
tion Ill. 


2. Corner transfer matrix renormalization group 

The corner transfer matrix renormalization group 
(CTMRG) was first introduced by Baxter [42]. It was 
further developed and applied to classical statistical sys¬ 
tems by Nishino and Okunishi [43, 44]. More recently, 
it was adapted to the contraction of tensor networks by 
Orus [33, 45]. CTMRG determines the ‘environment ten¬ 
sor’ Q of the four sites A, B, C, D as defined in Fig. 7. The 
locations of these four sites correspond to the locations 
of the ‘impurity sites’ in TERG. The relation between 
the environment tensor and the four-spin reduced den¬ 
sity matrix will be given below. 

Similarly to the TERG, one starts from and (see 
Eq. (12)) located at the corresponding sites of the tensor 
network with the exception of the four sites A, B, C, D 
(see Fig. 7). After the complete contraction of this tensor 
network, one obtains twelve tensors Ci, Tf, Tf, C 2 , 

T|, C 3 , Tg, Tg , C 4 , T 4 , T 4 shown on the right side of 
Fig. 7. They constitute the environment tensor. In order 
to determine them, the CTMRG algorithm successively 
contracts more and more tensors from the network into 
these twelve tensors (see Fig. 7). 




FIG. 7. Tensor network used in CTMRG (left) and environment 
tensor (right). Here, there are no tensors at the sites A, B, C, D. 

In order to prevent exponential growth of the virtual 
bond size a renormalization is performed at each step 
just like in TERG. However, the details of these renor- 
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malization steps are somewhat different. In CTMRG the 
twelve tensors are renormalized by left, up, right and 
down ‘moves’ defined and described in the following four 
steps. We describe left moves only as illustrated in Fig. 8 , 
the others are done analogously. The description of the 
steps follows Orus [33]: 

Step 1. Insertion: insert two sets (columns) of ten¬ 
sors as shown in Fig. 8 . (The insertion of two sets is 
only necessary because of translational symmetry break¬ 
ing discussed in the previous subsections.) 

Step 2. Absorption: absorb the first set of tensors 
into new tensors with increased vertical bond size: C( = 
CiT^, Tl' = r|T“, Tf = T^T\ C'^ = C^T^. (Here and 
in the following we omit the indices of the tensors, since 
they are easily reconstructed from the corresponding fig¬ 
ures.) 

Step 3. Renormalization: Insert two types of approx¬ 
imate isometries Z {Z^Z pe I) and W {W^W ~ I) as 
shown on the left side in Fig. 8 such that the vertical 
bond size of the tensors C[ , Tl', Tf, is truncated. The 
renormalized tensors are Ci = T| = ZT^'W^, 

T 4 = WT^'Z^, C 4 = ZC'^, and I is the identity matrix. 
One determines Z from an eigenvalue decomposition of 




FIG. 8. Left move of a CTMRG transformation. The tensors in 
the shaded box (right part of the figure) are inserted into the envi¬ 
ronment tensor (step 1). The inserted tensors are then combined 
with the left column of tensors (absorption, step 2) and renormal¬ 
ized (step 3) as illustrated in the center figure. The renormalized 
tensors shown in the left part of the figure replace the boxed part of 
the right part of the figure to form the renormalized environment 
tensor. 


the matrix C[Ci + C'^C'^ = ZDzZ^ and W from an 
eigenvalue decomposition of the matrix Q'lQi +Q 4 Q 4 = 
WDwW^ with Q[ = C[T^', Q '4 = In order to 

achieve the desired truncation (Z —Z, IT —>■ W) one 
keeps only the eigenvectors belonging to the Dc largest 
eigenvalues oi Dz and Dw, respectively. 

Step 4- Repeat steps 2 and 3 for the second set of 
inserted tensors. 

After the absorption and renormalization of the sec¬ 
ond set of tensors, one obtains the renormalized tensors 
Cl, r|, T 4 , (74 for the left column of tensors of the envi¬ 
ronment. A sequence of one left, down, right, and up 
moves constitutes one CTMRG transformation. Obvi¬ 


ously, this transformation resembles a coarse-graining of 
the tensor network. The moves described above are re¬ 
peated until convergence is achieved. We use the same 
convergence condition as for TERG with four-spin re¬ 
duced density matrix given by 


= Tr{gT^T^T^T% (19) 


in terms of the environment tensor Q and the unrenor¬ 
malized tensors T^, T^, T®, defined in Eq. (13). Of 
course, the reduced density matrix has to be normalized 
such that Trp = 1 . Alternatively, one may renormalize 
until for some n: ~ where A^ 

is the singular matrix of the corresponding corner tensor 
C^ [46]. 

In order to start up the recursive renormalization de¬ 
scribed above, all 12 tensors constituting the environ¬ 
ment are set to tensors and T^, respectively, and su¬ 
perfluous indices are traced out. 

A comparison of the two results Eq. (18) and (19) for 
the four-spin reduced density matrix may be instructive. 
In TERG the renormalized ‘impurity tensors’ T* con¬ 
tain all information about the tensor network, while in 
CTMRG one finally has to contract the unrenormalized 
impurity tensors into the renormalized ‘environment ten¬ 
sor’ Q in order to get the density matrix. The renormal¬ 
ization procedures used in both methods in order to pre¬ 
vent exponential growth of indices are somewhat differ¬ 
ent, however, it becomes obvious from the above descrip¬ 
tions that the two methods are in fact closely related. 


D. Translationally invariant matrix product states 

In this section we will briefly describe the methods we 
use for ID systems. Calculations in ID are numerically 
far less demanding than 2D calculations. However, it is 
instructive to compare different methods. 

The most efficient methods for ID calculations are 
variational methods. They are described in detail by 
Schollwdck [30]. Various imaginary time evolution al¬ 
gorithms have also been investigated for ID, notably the 
TEBD algorithm proposed by Vidal [36] . This algorithm 
motivated the 2D algorithm described in section HA. As 
discussed there, the TEBD method locally breaks trans¬ 
lational invariance. 

Here, however, we would like to discuss a method which 
maintains translational invariance exactly, i.e. we repre¬ 
sent a state j^k) of N spins by a matrix product state 
(MRS) with identical matrices at each lattice site 

1^)= X! Tr(A‘^i •...•A‘^^)lfTi,...,ajv). (20) 

(7l,...,(TAr 

The rank-3 tensors (A^ p,) have physical (spin) index 
cr of size two (since we consider spin- 1/2 systems only) 
and virtual dimensions of size m. Such an MPS was al¬ 
ready introduced in the seminal papers by Ostlund and 
Rommer [34, 35]; the PEPS introduced in Eq. (2) is its 























straightforward 2D generalization. We assume periodic 
boundary conditions. 

In order to implement imaginary time evolution with¬ 
out locally breaking translation invariance one requires 
a matrix product operator (MPO) representation of the 
time evolution operator exp(—riJ), 

exp(-TiJ) ^ Tr (r) • (r) • . ■. • 

CTl , . . . , (TiV , 

|cri,...,crAr)(cri,...,cr^|, (21) 

with the physical bonds a and a'. The trace is taken 
over virtual indices. The size of the (virtual) dimensions 
of the rank-4 tensors {W{[F ) depends on the details of the 
evolution operator under consideration and will be deter¬ 
mined for specific cases in the following. Details on MPO 
representations and their practical use may be found in 
Ref. [30]. Of course, in ID we also assume Hamiltoni¬ 
ans with nearest neighbor interactions only, and the same 
considerations about Trotter expansions as in the 2D case 
apply, i.e. the time evolution will proceed in small time 
steps At. 

Application of an MPO to an MPS will produce an 
MPS in terms of matrices A''^ with increased virtual bond 
dimension, 

KipW'p’)=T.^w'^pp'^ ( 22 ) 

(T' 

and in order to prevent an exponential growth of the MPS 
size at each step of imaginary time evolution, we need to 
truncate the size of the MPS at each evolution step. 

In order to do so one projects the MPS matrices A''^ 
to matrices of the same size as the original matrices 
As suggested in Ref. [47] one may use a projection op¬ 
erator which is similarly constructed as in DMRG using 
the density or transfer matrix, E = ® ■ The 

leading eigenvector P of A is rewritten as a square ma¬ 
trix, and from its singular value decomposition only the 
lowest m singular values are kept. This defines a projec¬ 
tion operator P which is used to project A"^ to a matrix 
A'^ = P^A"^P of the same dimensions as the original 
matrix A^. However, A^ corresponds to a later time of 
the system’s evolution. This renormalization procedure 
is illustrated in Fig. 9. After many imaginary-time steps 
and occasional reduction of the step size one reaches an 
approximate MPS representation of the ground state of 
the interacting spin system. 

For some time evolution operators, MPO representa¬ 
tions can be determined exactly. E.g. for the interaction 
of a spin with an internal or external field we use the 
identity {i = x,y,z) 

gKCTi _ cosh(/t)l-I-sinh(H;)tTi. (23) 

which can be proved using the properties of the Pauli 
matrices af ='I. Here 1 denotes the identity matrix. 



-d)- = 


FIG. 9. MPS renormalization 

The MPO representation for the evolution operator 
gKjZfc jg easily obtained in terms of 

the I X 1 X 2 X 2 tensors 

W = ^cosh(K)l -t-sinh(K)cri^ , (24) 

i.e. the size of virtual dimensions is 1. We have written 
the W tensor in terms of a variable k = rg, where g is 
the coupling strength of the field under consideration. 

For spin-spin interactions we need the identity 

gK(Ti®(Ti _ 0 1 -I- sinh(K)((Ti 0 (Ti). (25) 

With this relation one easily finds an MPO representa¬ 
tion of the evolution operator - 

in terms of the 2x2x2x2IF tensors, 

cosh(K)l a/ sinh(K) cosh(K)(Ti \ 

sinh(«;) cosh(K)(Ti sinh(«;)l J ’ 

(26) 

the size of the virtual dimension is 2. The latter relation 
was derived using a slightly different notation in Ref. [47] . 

The projection procedure for the reduction of the MPS 
size after each imaginary time step is the computationally 
most expensive part of the calculations to be performed. 
Therefore, it is desirable to streamline this step as much 
as possible. In fact it is desirable (from the computa¬ 
tional viewpoint) that the W tensors are real symmetric. 
This then leads to a symmetric transfer operator, which 
can be diagonalized rather efficiently. Suitable proce¬ 
dures for the symmetrization of W tensors are discussed 
in Ref. [47]. 

Our present realization of the translationally invariant 
MPS algorithm with real symmetric tensors W is appli¬ 
cable only for nonnegative values of k. In order to an 
have efficient algorithm for negative k one has to derive 
additional real W matrices from Eq. (25) (otherwise com¬ 
plex numbers appear intrinsically in W). Instead, we use 
the standard TEBD algorithm [25] in our calculations for 
parameter dependencies which correspond to negative k. 
The description of the TEBD algorithm is present in var¬ 
ious papers (e.g., [36]), thus we do not provide it in the 
present paper. 

For the translationally invariant MPS algorithm physi¬ 
cal quantities are calculated from the 2-spin reduced den- 






9 


sity matrix 

= Tl'iG ' (27) 

in terms of the environment matrix Q = {v l ® v/\^ 
and the unrenormalized ‘impurity matrices’ T’^ = 

® A'^. The environment matrix is determined from 
the leading left vl and right vr eigenvectors of the trans¬ 
fer matrix and its eigenvalue A. In this way 

we obtain results in the thermodyiramic limit (see [30]). 
Various results calculated from this density matrix are 
compared with 2D results in the next section. 

In the case of TEBD algorithm, a bipartition in the 
state representation is present, and the translationally in¬ 
variant ground state is represented by two MPSs {A, B}. 
The 2-spin reduced density matrix is calculated as 

=Tr(a-r;i"^-T^^"=). (28) 

The impurity matrices are here: = A'^ * <S> 

A'^, Tg'°' = 3"^'* ® B"^. The environment matrix is t/ = 
(wl®w_r)^/ a with A the leading eigenvalue and vl, vr the 
corresponding eigenvectors of a combined two-site trans¬ 
fer matrix (E., • (E., 

III. Entanglement measures and entanglement 
distribution: Numerical results and physical 
interpretation 

In this section we apply the formalism presented in 
the previous section to ID and 2D spin-I/2 systems: the 
quantum Ising model in a transverse magnetic field and 
the XXZ model. We calculate various entanglement mea¬ 
sures for these systems such as one-site entanglement en¬ 
tropy, one-tangle, concurrence of formation and negativ¬ 
ity. Furthermore, we determine bounds on the localizable 
entanglement in terms of the concurrence of assistance 
and maximal two-point correlation functions, local en¬ 
tanglement, and entanglement per bond. We compare 
these quantities and discuss their ability to identify crit¬ 
ical points and distinguish between different phases. En¬ 
tanglement per bond is presented only for 2D models, due 
to the fact that our translationally invariant MPS algo¬ 
rithm does not provide the MPS in canonical form [25]. 
The mentioned entanglement measures are briefly defined 
in Appendix IV. 

For all calculated quantities, the numerical results ob¬ 
tained from the TERG and CTMRG methods are nearly 
identical. Differences between both methods increase 
slightly in the critical region, and are strongly depen¬ 
dent on cutting parameters used in the renormalization 
procedures. A more complete analysis of such issues is 
under way. 

An interesting characteristic we analyze using the cal¬ 
culated entanglement measures is the monogamy of en¬ 
tanglement [9] or - more precisely - the entanglement 


distribution between different parties. Somewhat naively, 
entanglement monogamy may be expressed as follows: if 
two parties are maximally entangled they cannot be en¬ 
tangled at all with a third party. Expressions for the 
distribution of entanglement in the form of monogamy re¬ 
lations for multi-qubit systems, based on the concurrence 
of formation Cr and concurrence of assistance Ca have 
been obtained in Refs. [48, 49]. Thus, among the different 
entanglement measures we calculate in the present work, 
concurrence of formation and assistance are of primary 
interest. For the models studied in the present paper, 
a naive entanglement monogamy analysis was done in 
Ref. [15] using Monte-Carlo methods for the calculation 
of C'_F. Here, we provide a more comprehensive analysis 
based on monogamy relations for Cr and Ca- 

Entanglement monogamy relations for A^-qubit sys¬ 
tems are obtained from the Coffman-Kundu-Wootters 
(CKW) [9] inequality, 

- [^^‘’]aBi + [C'f]aB2+-■ • + 

(29) 

where [C'f]^^^ = [CR]{pABi) is the concurrence of the 
reduced density matrix pABi and [Cr]a\BiB 2 ...Bn-i ~ 
C(I'0 )a|BiB2...Sa,_i) the concurrence of the pure state \tp) 
as defined in Ref. [50] . A more general form of the CKW 
inequality was recently suggested in Ref. [51] . For N-qubit 
states, Bn-i obtained from the one-site 

reduced density matrix, and it is equal to the one-tangle 
d: [GF?AlB^B 2 ...B^_^ = 2(1-Trp^) = 4detpA = n [50]. 

In our analysis we use two main assumptions concern¬ 
ing the entanglement structure of the ground states of the 
models we study. The first is that only the nearest neigh¬ 
bor concurrences give major contributions to the sum of 
the right hand side of (29). The larger the separation 
between two particles the smaller is the concurrence be¬ 
tween them. The second assumption is a consequence of 
the translational symmetry of the ground states and as a 
consequence all nearest neighbor concurrences are equal. 

Taking into account these two features of the systems 
under consideration allows us to rewrite the inequality 
(29) for ID and 2D models. For ID systems one obtains 

+ (30) 

where [C'^'^] is the nearest neighbor concurrence of for¬ 
mation, and the quantity contains all other bipartite 
concurrences. Analogously in 2D one finds 

rr>4[ClP]l + Sr. (31) 

From these relations we obtain information about the en¬ 
tanglement distribution in the state. A more complete 
analysis of the entanglement distribution requires taking 
into account next-nearest neighbor two-party and longer 
ranged bipartite terms in the GKW-inequality. More¬ 
over, it is also possible to calculate three-party entangle¬ 
ment terms and look for their contribution to the entan¬ 
glement distribution. This is numerically easily feasible 
for ID systems, but is much harder for 2D models. 
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Dual to the CKW inequality one can derive the follow¬ 
ing relation wich involves the concurrence of assistance 
Ca on the right hand side [48, 49] 

- [^a\\bx + V^a]\b2.~^---'^[(^aTaBn-i- 

(32) 

Again we introduce the quantity 

+ (33) 

where contains the nearest neighbor terms and 

5a the longer-ranged bipartite concurrences. In 2D we 
have 

rF<^[CT]L+5T. (34) 

A. Quantum Ising model in a transverse field 

The spin-i Ising model in a transverse magnetic field 
h is given by the Hamiltonian 

Rising = J ^ ^ ^ ^ ( 35 ) 

<bi> ® 

where the Ui (i = 1,2,3) are the standard Pauli spin 
operators. This model is symmetric (spin-flip sym¬ 
metric) . 

The sign of the coupling constant J determines the 
type of the interaction between spins: anti-ferromagnetic 
for J > 0 and ferromagnetic for J < 0. The calculated 
physical quantities are symmetric with respect to J = 0. 
Quantities like the magnetization {J < 0) are mapped 
to their staggered counterpart (J > 0). In the present 
paper we choose the energy scale by setting J = — 1. 

In ID this model can be solved analytically using a 
Jordan-Wigner transformation [52]. It is well known that 
at the critical points h = ±1 this model shows quan¬ 
tum phase transitions separating a magnetically ordered 
phase (—1 < h < I) from paramagnetic phases {h < —1 
and h > 1). In the ordered phases the Z 2 symmetry is 
spontaneously broken. At the critical points and in the 
thermodynamic limit the ground state energy per site is 
given by Eq = — 4 / 7 r. 

The 2D quantum Ising model cannot be solved ana¬ 
lytically, and various methods are applied to solve it nu¬ 
merically, notably rather resource-intensive Monte-Carlo 
(MC) methods. Such calculations find a transition be¬ 
tween a ferromagnetic and a paramagnetic phase at a 
critical point her = 3.044 [53]. The tensor network im¬ 
plementation we use here produces numerical results sig¬ 
nificantly faster than MC calculations, however, with less 
precision: our implementation determines a critical point 
at her ~ 3.25, which is determined from a singular point 
of the second derivative of the ground state energy as a 
function of h. Of course, significantly more precise results 
could be obtained with more elaborate tensor network 
implementations and larger bond sizes [23]. However, it 


is our goal to investigate correlations and entanglement 
properties using small numerical cost. In practice, we 
study positive h only and obtain results for negative h by 
a reflection at h = 0. In order to compare numerical re¬ 
sults for one- and two-dimensional systems we rescale the 
magnetic field dependence h/hcr such that phase transi¬ 
tions always occur at h/hcr = 1- 

For h ^ he, the system (both in ID and 2D) is a 
classical Ising model with a doubly degenerate ground 
state (in the thermodynamic limit). In experimental sit¬ 
uations this degeneracy is broken and this is done intrin¬ 
sically in our MPS and PEPS implementations as well. 
For h ^ her the magnetic field dominates and the ground 
state corresponds to free spins oriented according to the 
magnetic field. 

With Fig. 10 we start the presentation of the numerical 
results and show the magnetizations and as a 
function of h. We see that m^ih/hcr) in ID increases 
slower from zero magnetic field towards the critical point 
and has lower value at the critical point then the ttIt in 
the 2D. 

At this stage we do not quantitatively extract critical 
exponents as this would require more precise and time- 
consuming calculations close to the critical points. How¬ 
ever, qualitatively the critical properties are in agreement 
with expectations. 



h/hcr 


FIG. 10. (colour online) Magnetizations nix and niz as a function 
of the magnetic field h/hcr for the ID and 2D quantum Ising mod¬ 
els. Parameters tor the ID MPS calculation: m = 20. Parameters 
for the 2D TERG calculation: D = 4, Dc = 20. 

In Fig. 11 we show the entanglement measures calcu¬ 
lated from the single-spin density matrix: one-site en¬ 
tanglement entropy Si and one-tangle ti for the one- 
and two-dimensional Ising models. It is clearly seen that 
all these measures nicely peak in cusps at the critical 
point. All 2D results are multiplied by a factor of 2 for 
the convenience. 

In Figs. 12 and 13 we show entanglement measures cal¬ 
culated from the two-spin density matrix as a function of 
the magnetic field: the concurrence of formation Cp and 
negativity N. Our TN results are in a very good agree¬ 
ment with the Monte Carlo results by Syljuasen [15]. Of 
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h/hcr 


FIG. 11. (colour online) One-site entanglement entropy Si and 
one-tangle ri as a function of the magnetic field h/hcr for the ID 
and 2D quantum Ising models. 2D results are multiplied by a factor 
of 2. Parameters for the ID MPS calculation: m = 20. Parameters 
for the 2D TERG calculation: D = A, Dc = 20. 


course, calculations close to the critical point in 2D are 
difficult both for TN and MC methods, but clearly the 
cusp at the critical point can be better resolved with the 
TN method used here. Close to the critical point the 
MC results of Ref. [15] are very noisy. The concurrence 
of formation for ID does not peak at the critical point, 
but shows an inflection. This is in agreement with the 
analytical results presented in Ref. [8]. 

The negativity shows similar characteristics as the con¬ 
currence of formation both in ID and 2D. Concurrence 
of formation and negativity reach their maximum at the 
same value for the magnetic held. Negativities for ID 
and 2D geometries both satisfy the concurrence bounds 
V'(l - Cf)'^ + C'j, - (I-Cf) < N <Cf [54]. 



h/hcr 


FIG. 12. (colour online) Concurrence of formation Cf and nega¬ 
tivity Af as a function of the magnetic field h/hcr for ID and 2D 
quantum Ising model. 2D results are multiplied by a factor of 2. 
Parameters for the ID MPS calculation: m = 20. Parameters for 
the 2D TERG calculation: 0 = 4, Dc = 20. 

The local entanglement 5'ioc shown in Fig. 13 is the 
simplest form of a block-block entanglement, the entan¬ 


glement between two neighboring spins and their envi¬ 
ronment. We see that both in ID and 2D this measure 
has a peak with a cusp at the critical point. Not sur¬ 
prisingly, S'loc’s absolute value at the critical point is the 
largest among other entanglement measures we calculate 
from the two-site reduced density matrix. This is due to 
the fact that S'loc correspond to entanglement between 
two neighbor spins as one party with all other spins as 
another party in contrast to the entanglement between 
just two neighbor spins in the case of Cf, N, Ca- Sim¬ 
ilar to the one-site entanglement entropy, S'loc is small 
in the paramagnetic phase {h < her), increases sharply 
close to the critical point and then decreases slowly in 
the ferromagnetic phase {h > her)- 

Fig. 13 also demonstrates that the bipartite entangle¬ 
ment per bond in 2D identihes the critical point hav¬ 
ing a peak with a cusp there. This measure exemplihes 
one useful advantage of the translationally invariant TN 
methods: the possibility to extract information about the 
state right from the TN representation, i.e. one does not 
need to calculate expectation values at potentially high 
numerical cost. 



h/hcr 


FIG. 13. (colour online) Comparison of local entanglement Sio^ 
dependence on magnetic field h for the ID and 2D quantum Ising 
model. Entanglement per bond Sps dependence for the 2D Ising 
model. Results are renormalized to the h/hcr dependence. Results 
for the 2D model are multiplied by a factor of 2. Parameters for 
the ID MPS calculation: m = 20. Parameters for the 2D TERG 
calculation: D = 4, Dc = 20. 

In Fig. 14 we compare the upper bound (concurrence of 
assistance Ca) and lower bound (maximal two-site corre¬ 
lation function Qmax) of the localizable entanglement as 
a function of the magnetic held. Our results show cusps 
at the critical point and are in a good agreement with 
those obtained using other methods [15, 55]. Note, that 
our MPS and PEPS implementation intrinsically break 
the Z 2 symmetry, thus leading to product states for small 
and large magnetic helds. 

All entanglement measures discussed above are able to 
identify the critical point of the system both in ID and 
2D. The fastest and easiest way to identify the critical 
point is obtained from the entanglement per bond. This 
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FIG. 14. (colour online) Upper and lower bounds of the localizable 
entanglement as a function of magnetic field h/hcr for the ID and 
2D quantum Ising models; 2D results are multiplied by a factor of 
2. The shaded areas between Ca and Qmax for ID and 2D results 
correspond to possible values of the localizable entanglement. Pa¬ 
rameters for the ID MPS calculation: m = 20. Parameters for the 
2D TERG calculation: D = A, Dc = 20. 


measure explicitly requires a tensor network representa¬ 
tion and cannot be obtained using other methods. As 
expected, all entanglement measures approach zero for 
small and large transverse magnetic fields, which indi¬ 
cates product states for these limits. 

In Fig. 15 we show the concurrence of assistance 
2 the one-tangle , and the concurrence of 

formation 2 Comparing t™ and 2 we 

see that the CKW inequality (29) is fulfilled and that the 
nearest-neighbor two-particle entanglement corresponds 
to only about 25% of the bipartite entanglement in the 
critical region. At the same time, outside of the critical 
region 2 \C]P'\^^ nearly exhausts the CKW inequality. 
This behaviour quantitatively confirms that the phase 
transition is characterized by the presence of long-range 
entanglement. Comparing and (2 we con¬ 

clude that already the nearest neighbor entanglement 
contributions are larger than the lower bound of how 
much entanglement can be created by assistance. 

Fig. 16 displays the entanglement monogamy analy¬ 
sis for the 2D quantum Ising model. Here, we compare 
and 4[C|P]^^.. The CKW inequality is 
fulfilled, and the nearest-neighbor entanglement in the 
critical region corresponds to about 50% total bipartite 
entanglement. In comparison to the ID result, we ob¬ 
serve that the 2D nearest-neighbor entanglement con¬ 
tains more of the total bipartite entanglement, which 
can be explained by the presence of a larger number 
of nearest neighbors of each site. Again, similar to the 
ID case, 4 nearly exhausts the CKW inequality 

outside of the critical region. Again, comparing and 
4 [C\% . we see that also in 2D the nearest-neighbor en¬ 
tanglement terms in general are already larger than the 


FIG. 15. (colour online) Entanglement monogamy analysis for 
the ID quantum Ising model: Comparison of the concurrence of 
formation Cf, the concurrence of assistance Ca the 1-tangle 
ri. For details see the discussion in the main text. Parameters for 
the ID MPS calculation: m = 20. 


lower bound on how much entanglement can be cre¬ 
ated by assistance. 



h 

FIG. 16. (colour online) Entanglement monogamy analysis for 
the 2D quantum Ising model: Comparison of the concurrence of 
formation Cp, the concurrence of assistance Ca ^-ud the 1-tangle 
ri. For details see the discussion in the main text. Parameters for 
the 2D TERG calculation: D = A, Dc = 20. Critical value of the 
magnetic field is her ~ 3.28. 


B. XXZ model 

Next we study the spin-i XXZ (anisotropic Heisen¬ 
berg) model, 

jjXXZ = I (g) (jJ -1- erf 0 erj -I- A (jf 0 cr|} (36) 

(rj) 

as a function of the anisotropy parameter A. The Hamil¬ 
tonian of this model is U(I) symmetric (corresponing to 
an invariance under a U(l) rotation about the spin ^ 
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axis) as well as Zi symmetric (corresponding to an in¬ 
variance under a tt rotation about the spin x or y axis). 
It is SU(2) symmetric at the Heisenberg point A = 1. 
The ground state of the XXZ model in different phases 
preserves these symmetries or not depending of the space 
dimension. [56]. The Z 2 symmetry implies that (erf) = 0 
and {(jfo-p = (efo-j) = 0. The U(l) symmetry implies 
that «) = (af) = 0 , (afap = (afaj), (afaj) = 0 . 

The XXZ model has a richer phase structure than 
the Ising model: The ID XXZ model shows three 
phases [57, 58]. For A > 1 the system is in a gapped 
antiferromagnetic phase (in particular, it corresponds to 
a classical Ising anti-ferromagnet for large positive A). 
At A = 1 there is a critical point, where an infinite- 
order Kosterlitz-Thouless quantum phase transition oc¬ 
curs from the anti-ferromagnetic phase to the XY phase. 
The system is equivalent here to the spin-^ Heisenberg 
anti-ferromagnet with a gapless ground state. In the XY 
phase (JAj < 1 ) the system is gapless and the correlation 
functions decay polynomially. At A = — 1 the system un¬ 
dergoes a first-order quantum phase transition to a fer¬ 
romagnetic gapped phase for A < — 1. For large negative 
A the system resembles an Ising ferromagnet. 

In the thermodynamic limit spontaneous Z 2 symme¬ 
try breaking occurs in the ferromagnetic (A < —1) and 
antiferromagnetic phases (A > 1), but Z 2 symmetry is 
preserved in the XY phase. The continuous U(l) symme¬ 
try remains unbroken in all three phases of the ID XXZ 
model. 

The two-dimensional XXZ model shows three different 
phases, as well [59-61]: an antiferromagnetic phase for 
A > 1, an XY phase for jA] < 1 and a ferromagnetic 
phase for A < —1. It undergoes a second-order phase 
transition at A = I [62] and a first-order phase transition 
at A = —I [63]. Just as in ID, the Z 2 symmetry is 
spontaneously broken in the ferromagnetic (A < — 1) and 
antiferromagnetic phases (A > 1 ), and remains unbroken 
in the XY phase. However, unlike in the ID case, the 
continuous U(l) symmetry is not preserved in the XY 
phase of the 2D XXZ model [56]. 

The one-dimensional XXZ model has been studied ex¬ 
tensively using the Bethe Ansatz [64-67] . At the Heisen¬ 
berg point (A = 1) the ground state energy is found to 
be Eq = —4 log 2 -1-1. In order to demonstrate the de¬ 
pendence of our numerical results on the MPS virtual 
dimension, in Table I we compare the ground state en¬ 
ergy calculated for different m with the analytical value. 
In the following we will present results calculated with 
TO = 20 which provides accurate results at moderate nu¬ 
merical costs. For the A < 0 region we will use the 
TEBD imaginary time evolution method instead of the 
translationally invariant MPS method, because our MPS 
implementation is not optimal for calculations in this re¬ 
gion by construction. Tests with both algorithms in the 
region A > 0 show that the results agree to a high pre¬ 
cision. 

The Heisenberg point (A = 1) for the 2D XXZ model 
was intensively studied in Ref. [ 68 , 69]. The best quan- 


m 

E 

AE 

10 

-1.77202 

1.0 10”^ 

15 

-1.77237 

3.8 10-'‘ 

20 

-1.77247 

2.1 10“^ 

25 

-1.77253 

1.1 lo-'^ 

30 

-1.77254 

8.5 10“® 


BA -1.77259 

TABLE I. Ground state energy Eq of the ID XXX model com¬ 
pared to the Bethe Ansatz (BA) result as a function of MPS virtual 
dimension m. The relative difference is AE = {Eq — £^ ba )/-^ ba - 

D Da E AE 

'2 20 -1.318 2.2 10”^ 

3 20 -1.327 1.3 10"^ 

4 32 -1.333 7.0 10"^ 

5 64 -1.338 2.0 10"® 

QMC TM) 

TABLE 11. Ground state energy Eq of the 2D XXX model com¬ 
pared to QMG results as a function of PEPS virtual dimension 
D and TERG cutting parameter D^. The relative difference is 
AE = (Eq — Eqmc)/15qMC' 

turn Monte Carlo result for the ground state energy is 
EQMC = —1.340 [69]. In Table H we compare our nu¬ 
merical results for this energy for different D and Dc- 
In the following we will present results calculated using 
D = A and Dc = 20. With this choice we obtain good 
results at reasonable numerical cost. 

In Fig. 17 we show the ground state energy per site 
as a function of the asymmetry parameter A for ID and 
2D. From the figure we see that in the parameter region 
A < — 1 the ground state energy is linearly dependent 
on A: Ao(A) = A for both the ID and 2D XXZ models. 
At A = 1 the ground energy shows a kink for 2D, but 
for ID is continuous and infinite-order differentiable as is 
known from analytical analysis [64-66]. 



A 

FIG. 17. (color online) Comparison of ground state energy as a 
function of A for the ID and 2D XXZ models. Parameters for 
the ID MPS calculation: m = 20. Parameters for the 2D TERG 
calculation: D = 4, Dc = 20. 

Fig. 18 shows various magnetizations as a function of 
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the asymmetry parameter A. Non-zero ruz magnetiza¬ 
tion in the ferromagnetic phase (A < —1) and non-zero 
staggered magnetization mf" in the anti-ferromagnetic 
phase (A > 1) both in ID and 2D models confirm the 
Z 2 symmetry breaking in these phases. In the XY phase, 
Fig. 18 shows a U(l) symmetry breaking not only for 
the 2D model, as expected, but also for the ID model. 
This clearly shows a deficiency of the numerical method 
used. This U(l) symmetry breaking is strongly depen¬ 
dent on the chosen m, and gets smaller with increasing m, 
however, it appears that one needs to use a code which 
implements the U(l) symmetry of the states from the 
outset in order to get more precise results. We will do 
this in a future paper. This unphysical breaking of the 
U(l) symmetry will be seen as well later in various cal¬ 
culated entanglement measures. Note that the transla- 
tionally invariant MPS algorithm nicely obtains the anti¬ 
ferromagnetic phase, despite the fact that it uses equal 
tensors at each site. The magnetizations in one dimen¬ 
sion are weaker than their two-dimensional counterparts. 




A 


FIG. 18. (colour online) Comparison of magnetization rUz and 
its staggered counterpart mj* as a function of A for the ID and 
2D XXZ models. U(l) symmetry breaking for |A| < 1 also gives 
nonzero staggered magnetization m®*. Parameters for the ID MPS 
calculation: m = 20. Parameters for the 2D TERG calculation: 
D = 4, Dc = 20. 

In Fig. 19 we show the one-site entanglement measures: 
one-site entanglement entropy Si and one-tangle ri for 
the one- and two-dimensional XXZ models. All these 
measures peak in cusps at the critical point A = 1 and 
are zero for the A < — 1. At the Heisenberg point in 


the ID model the ground state is SU(2) symmetric and 
the one-site measures Si and ti approach their maximal 
possible values. Theoretically it is expected that these 
quantities equal to 1 throughout the XY phase in ID, 
but due to the U(l) symmetry breaking introduced in the 
algorithms as discussed above these quantities decrease 
while approaching the A = — 1 critical point. Again in 
ID we would obtain better results for larger m or by 
using a code which respects the U(l) symmetry from the 
outset. 



FIG. 19. (colour online) Comparison of one-site entanglement 
entropy 5i and one-tangle ri dependence on A for the ID and 2D 
XXZ models. Results for 2D model are multiplied by a factor of 2. 
Parameters for the ID MPS calculation: m = 20. Parameters for 
the 2D TERG calculation: 0 = 4, Dc = 20. 

In Fig. 20 we show the two-site entanglement measures: 
concurrence of formation Cp and negativity N for the 
one- and two-dimensional XXZ models. The concurrence 
of formation for the ID and 2D XXZ models was studied 
in Refs. [15, 16, 56], and our results are in a very good 
agreement. The figure nicely shows that Cp and N in one 
and two dimensions have maxima exactly at the critical 
point A = 1. It is known that Cp is related to the ground 
state energy [70]. We see that similarly to the ground 
state energy in Fig. 17, Cp and N show a maximum 
at the critical point. Our results correspond to the fact 
discussed in Ref. [61] that the ground state energy of the 
XXZ model in two and three dimensions shows a cusp at 
the transition point, thus leading to a cusp in concurrence 
of formation. The ID Cp and N just have maxima at 
the critical point A = 1 without cusps. 

Negativity for the XXZ model was previously stud¬ 
ied for a two-qubit chain [71] and for infinite tree ten¬ 
sor network states [72]. Our results extend such stud¬ 
ies to infinite chains and infinite square-lattice systems. 
Again, negativities for ID and 2D geometries both satisfy 
the concurrence bounds a/(1 — Cp)^ -I- Cp — (1 — Cp) < 
A^ < Cp [54]. Similarly to the quantum Ising model, 
we find that negativities for the ID and 2D XXZ models 
have a similar behavior as the concurrence of formation 
in ID and 2D. The 2D negativity peaks in a cusp and the 
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ID negativity just shows maximum at the critical point 
A = 1. 



A 


FIG. 20. (colour online) Concurrence of formation Cp and nega¬ 
tivity Af as a function of A for the ID and 2D XXZ models. Results 
for the 2D model are multiplied by a factor of 2. Parameters for 
the ID MPS calculation: m = 20. Parameters for the 2D TERG 
calculation: D = 4, Dc = 20. 

In Fig. 21 we present the local entanglement 5'ioc for 
the one- and two-dimensional XXZ models. The entan¬ 
glement per bond SpB for the 2D XXZ model is also 
shown. Local entanglement for the 2D XXZ model was 
studied in [17]. However, S'loc requires comment: for 
A ^ 1 we observe that the local entanglement reported 
here approaches zero while in Ref. [73] it approaches 1. 
The reason for this difference is the fact that the ground 
states we consider here has broken Z 2 symmetry, while 
the authors of Ref. [73] assume that the ground state is 
Z 2 symmetric. We see that local entanglement for ID 
and 2D has similar behaviour as the one-site entangle¬ 
ment entropy ^i. In both ID and 2D S'loc vanishes at 
A = — 1 and peaks in a cusp at A = 1. 

Entanglement per bond for the 2D XXZ model was 
analyzed in [18], but the authors discuss the SpB depen¬ 
dence on an external magnetic field with some fixed A. In 
our studies we have no external magnetic field and vary 
the anisotropy parameter A. Similarly to the quantum 
Ising model, SpB shows its ability to determine critical 
points by vanishing at A = —1 and having a peak with 
a cusp at A = 1. 

In Fig. 22 we show the upper bound (concurrence of 
assistance Ca) and lower lower (maximal two-site cor¬ 
relation function Qmax) on the localizable entanglement 
for the one- and two-dimensional XXZ models. These 
bounds in the two-dimensional case were also studied in 
Ref. [15]. 

We observe that for |A| < 1 the concurrence of assis¬ 
tance Ca and two-point correlation function Qmax de¬ 
crease for smaller A while in Ref. [15] Ca = 1 through¬ 
out the XY phase and Qmax does not drop to zero at one 
of the critical points. This difference in our results and 
results from [15] can be explained as following. It was 
shown in Ref. [56] that concurrence of formation Cp is 



A 

FIG. 21. (colour online) Local entanglement Sioc and entangle¬ 
ment per bond SpB as a function of A for ID and 2D XXZ model. 
Results for 2D model are multiplied by a factor of 2. Parameters 
for MPS calculation: m = 20. Parameters for TERG calculation: 
D = 4, Dc = 20. 


unaffected by spontaneous symmetry breaking (namely, 
U(l) symmetry breaking) for the zero-field XXZ-model. 
Let us consider also the concurrence of assistance Ca- 
The formula for Ca for maintained U(l) symmetry and 
broken Z 2 symmetry was introduced in [15]: 



Following the ideas from Ref. [56] for deriving the expres¬ 
sion for Cp for broken U(l) symmetry and maintained 
Z 2 symmetry, we find that Ca is in this case 

Ca = \ (^(l + (afaj))2-4(af)2 + l - «aj)) . 

(38) 

Obviously, Ca (unlike Cp) is affected by U(l) symme¬ 
try breaking. 

When U(l) and Z 2 symmetries are obeyed, Ca = 1- 
This is theoretically predicted, e.g., for the Heisenberg 
point A = 1. We see from our ID results that indeed 
Ca(A = 1) Ri 1. At the same time our 2D results for 
Ca for a = 1 do not reach the value Ca = 1- This can 
be explained by the fact that it is numerically hard to 
converge to the point where both (uf) and (erf) are zero, 
thus giving Ca = 1 from both equations (37) and (38). 

The discrepancy of our result for Qmax a-nd the corre¬ 
sponding result from Ref. [15] can be explained by U(l) 
symmetry breaking, resulting in a nonzero (erf). While 
(erf) increases, the function = (erferf) - (erf)(erf) 
(which is larger than and in the XY phase) de¬ 
creases to zero. 

Thus, we see that all entanglement measures discussed 
above are zero for A < — 1 and also approach zero for 
large positive A, indicating a product state in this limit. 
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FIG. 22. (colour online) Bounds on the localizable entanglement 
as a function of A for the ID and 2D XXZ models. Results for 
the 2D model are multiplied by a factor of 2. The shaded areas 
between Ca ^nd Qmax for ID and 2D results correspond to possible 
values of the localizable entanglement. Parameters for the ID MPS 
calculation: m = 20. Parameters for the 2D TERG calculation: 
D = A, Dc = 20. 


For the monogamy analysis in the Fig. 23 we repre¬ 
sent nearest neighbor entanglement, given by concur- 
rence of assistance (2 one-tangle and near¬ 

est neighbor entanglement, given by concurrence of for¬ 
mation (2 \pW\^nri)- comparing and (2 
we see that CKW inequality is fulfilled and the nearest 
neighbor two-particle entanglement corresponds to about 
1/3 fraction of the entanglement in the critical region 
around A = 1 critical point. For large A >> 1 the near¬ 
est neighbor two-particle entanglement approaches . 

By comparing and 2 we see that nearest 

neighbor entanglement in general is larger than the lower 
bound on how much entanglement can be created by 
assistance. Only in the XY phase in the region close 
to A = — 1 the nearest neighbor entanglement does not 
exceed the This feature is shown in the inset in the 
Fig. 23. 

Fig. 24 shows the entanglement monogamy analysis for 
the 2D XXZ model. In this case we compare (4 [0^^]^^.), 

and (4 The CKW inequality is fulfilled. 

Nearest neighbor two-particle entanglement is less then 
1/3 fraction of the entanglement in the entanglement 
distribution in the critical region around A = 1 criti¬ 
cal point. Similar to ID case, for high A ^ 1 nearest 
neighbor two-particle entanglement approaches 

By comparing and 4 we see that nearest 

neighbor entanglement in general is larger than the lower 
bound on how much entanglement can be created by 
assistance. And again, only in the XY phase in the region 
close to A = — 1 the nearest neighbor entanglement does 
not exceed ry, which is shown in the inset in the Fig. 24. 
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FIG. 23. (colour online) Entanglement monogamy analysis for the 
ID XXZ model: Comparison of the concurrence of formation Cp, 
the concurrence of assistance Ca and the 1-tangle ri. For details 
see the discussion in the main text. Parameters for the ID MPS 
calculation: m = 20. 



A 

FIG. 24. (colour online) Entanglement monogamy analysis for the 
2D XXZ model: Comparison of the concurrence of formation Cp, 
the concurrence of assistance Ca and the 1-tangle ti. For details 
see the discussion in the main text. Parameters for the 2D TERG 
calculation: D = 4, Dc = 20. 


IV. Conclusions 

We have investigated entanglement properties of infi¬ 
nite ID and 2D spin-1/2 systems using tensor network 
methods: the Ising model in transverse field and the 
XXZ model. Specifically we used a translationally in¬ 
variant MPS method in ID and TERG and CTMRG in 
2D in order to calculate the ground state of those mod¬ 
els. Different entanglement measures, such as one-site 
entanglement entropy and one-tangle, concurrence of for¬ 
mation and negativity, bounds on localizable entangle¬ 
ment (concurrence of assistance and two-point correla¬ 
tion function), local entanglement and entanglement per 
bond were calculated. 

Many of our results are in good agreement with those 
obtained using other numerical methods. This agreement 
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underlines that such tensor network methods are power¬ 
ful tools for the investigation of quantum models. The 
translationally invariant MPS algorithm and the TEBD 
algorithm lead to an U(l) symmetry breaking in the 
XY phase for the ID XXZ model and, therefore, our 
results are at variance with those assuming U(l) symme¬ 
try [57, 58]. 

Our results confirm the observation [18] that the bipar¬ 
tite entanglement per bond can successfully determine 
critical points. This measure is unique to tensor network 
methods. 

We made an entanglement monogamy analysis: The 
Coffman-Kundu-Wootters inequality is fulfilled in both 
models we studied, and the obtained entanglement dis¬ 
tribution indicates the presence of a relatively large frac¬ 
tion of long-range entanglement in the critical region for 
both Ising and XXZ models in both ID and 2D. 

Our work may be extended into several directions: In 
order to more deeply analyze the numerical possibilities 
of the translationally invariant MPS algorithm one needs 
to implement it efficiently for negative parameter values, 
that is negative k from Section II. Furthermore, it is de¬ 
sirable to have codes where the symmetries of the ground 
states can be prescribed from the outset. Such work is 
under way. 

Moreover, for more complete entanglement characteri¬ 
zation of the models it is important to take into account 
other entanglement measures and characteristics, such as 
fidelity [74], global entanglement [18], and entanglement 
spectrum [ 11 ]. Another promising direction is the anal¬ 
ysis of the complementarity of the entanglement [75] in 
many-body systems. And, of course, it would be inter¬ 
esting to extend our studies to higher spins. 

Appendix: Entanglement measnres 

In this appendix we briefly review well known defini¬ 
tions for various bipartite entanglement measures. 

The first two are the one-site entanglement entropy 
Si and one-tangle ri, which are obtained directly from 
the single-site reduced density matrix. The entanglement 
entropy [76] for bipartite pure states |' 0 i 2 ) is the von Neu¬ 
mann entropy of the reduced density matrix 

5(1^12)) =5(pi)=5(p2), (39) 

with the reduced density matrices pi = Tr 2 (pi 2 ) and 
P 2 = Tri{pi 2 ); pi 2 = I^i 2 )(^i 2 | and Tr^ indicates a trace 
over the subsystem i. The von Neumann entropy S' of a 
density matrix p is calculated from its eigenvalues [76] 
A.: 

S{p) =-plog2P= (40) 

i 

In the main text we use Si = S{pi). The one-tangle [9] 
is also calculated from one-site reduced density matrix: 

ri(pi) = 4detpi. (41) 


The von Neumann entropy is connected to the one-tangle 
through the relation [ 6 ] 

+ (42) 

where h{x) = —x log 2 a; — (1 — a;) log 2 (l — x) denotes the 
binary entropy function. 

Next we mention measures obtained from the two-site 
reduced density matrix pi 2 . A simple measure of bipar¬ 
tite entanglement in a mixed state is the entanglement 
of formation, Ep [77]. It counts the minimum number of 
maximally entangled states (Bell states) needed to con¬ 
struct a given state using only local operations and clas¬ 
sical communication (LOCC) (for details see [76, 77]). 
The entanglement of formation can be calculated from 
the concurrence of formation Cp [78, 79]: 

where h{x) denotes the binary entropy function. The 
coneurrence of formation [79] is an entanglement mea¬ 
sure for mixed states of two qubits, defined as 

Cp{p) = inax(0, Ai - A 2 - A 3 - A 4 ), (44) 

where Ai, A 2 , A 3 , A 4 are the eigenvalues in decreasing or¬ 
der of the Hermitian matrix 



with pi 2 = {ay ® ay)pl 2 {ay 0 ay). Here Pi 2 is the com¬ 
plex conjugate of the two-site density matrix pi 2 . Alter¬ 
natively, Xi are the square roots of the singular values 
of the non-Hermitian matrix pi 2 pi 2 - The concurrence is 
zero for a product state and one for a maximally entan¬ 
gled state. 

Another type of concurrence, the concurrence of assis¬ 
tance Ca, was introduced in connection with the entan¬ 
glement of assistance Ea [80]. Ca is obtained from [81]: 

Ca = Ai -I- A 2 -|- A 3 -I- A 4 . (46) 

The entanglement of assistance measures the maximal 
bipartite entanglement which be obtained while doing 
measurements on the rest of the spins. The idea of en¬ 
tanglement of assistance originates from the analysis of 
tripartite systems, described by a state By vary¬ 

ing the measurement on party 3, the ‘helper’ 3 is able 
to influence the mixed state of parties I and 2 [82]. In 
order to use Ea in practice one must be able to perform 
a maximization over all different measurement strategies, 
thus this measure is difficult to calculate. However, there 
exist easily calculable bounds on Ea'- upper bounds on 
Ea are the entropic bound, the fidelity bound, and con¬ 
currence bound Ca [80] . The latter is used in the present 
paper. 
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The localizable entanglement [55] is defined as the 
maximal amount of entanglement that can be localized 
(on average) between two spins while doing only local 
measurements on the rest of the spins in the environ¬ 
ment. El cannot be obtained from the reduced density 
matrix alone, thus it is able to describe characteristics 
of the wave function that are not captured by two-point 
correlation functions, e.g. exotic phases like topological 
orders. The calculation oi El is not a trivial task since 
one needs to optimize over all possible local measurement 
strategies, nevertheless it is possible to obtain bounds on 
El using only two-point correlation functions [55]. 

The upper bound for El is the concurrence of assis¬ 
tance Ca-, and the lower bound is obtained from the max¬ 
imal two-point correlation function, 

max(lQ?2"l, IQfll, \Qr2\) <El< Ca, (47) 

where (5i2^(|V')(V'l) = (V’k“ 0 crflV’) - (V’K 0 

l 2 |'*/’)(V^|li ® ll^) '^oc are the Pauli spin matrices. 

The negativity [83] is an ‘easy-to-compute’ measure de¬ 
fined as 

A7(pi2) = (48) 

where ^^"2 partially transposed density matrix pi 2 

with respect to subsystem 1. And IIP 12 II 1 = Try^P 12 P 12 
is the trace norm. IIP 12 II 1 is calculated as a sum of the 
singular values of pi 2 . A measure closely related to the 


negativity is the logarithmic negativity [84], 

En{p) =\og2{\\p^^\\i). (49) 

A simple form of bipartite entanglement is the entan¬ 
glement between two neighboring spins and the other 
spins of the system. This measure is called local en¬ 
tanglement [17]. The two-site local entanglement S'loc 
is obtained by tracing out all spin degrees of freedom of 
the system except the two nearest-neighbour spins and 
then calculating the von Neumann entropy of the result¬ 
ing reduced density matrix pi 2 , 

Sloe = S{pi2). (50) 

Another entanglement measure, which can be used if 
we have available a tensor network representation of the 
state in conventional form, is the bipartite entanglement 
per bond S'pb [18]. It is obtained from the bond vec¬ 
tors [25] (see also section 2) connecting two neighboring 
sites. The bond vectors contain essential entanglement 
information of the system. The entanglement per bond 
S'pb is given by 

SpB = ^A?log2A^ (51) 

i 

where the components of the bond vectors are normalized 
such that ^ ■ Af = 1. 

There are other entanglement measures like fi¬ 
delity [74], global entanglement [18], entanglement spec¬ 
trum [11] with Schmidt gap, which also can be used to an¬ 
alyze entanglement in many-body systems. These mea¬ 
sures are not considered in the present text. 


[1] S. Sachdev, Quantum Phase Transitions. 2nd Edi¬ 
tion (Cambridge University Press, Cambridge, England, 
2011 ). 

[2] X.-G. Wen, Quantum Field Theory of Many-Body Sys¬ 
tems (Oxford University Press, New York, USA, 2004). 

[3] H. T. Diep, Frustrated spin systems. 2nd Edition. (World 
Scientific, Singapore, 2013). 

[4] T. Sentliil, A. Visliwanath, L. Balents, S. Sachdev, and 
M. P. A. Fisher, Science 303, 1490 (2004). 

[5] T. Sentliil, L. Balents, S. Sachdev, A. Vishwanath, and 
M. P. A. Fisher, Phys. Rev. B 70, 144407 (2004). 

[6] L. Amico, R. Fazio, A. Osterloh, and V. Vedral, Rev. 
Mod. Phys. 80, 517 (2008). 

[7] T. J. Osborne and M. A. Nielsen, Phys. Rev. A 66, 
032110 (2002). 

[8] A. Osterloh, L. Amico, G. Falci, and R. Fazio, Nature 
416, 608 (2002). 

[9] V. Coffman, J. Kundu, and W. K. Wootters, Phys. Rev. 
A 61, 052306 (2000). 

[10] R. Horodecki, P. Horodecki, M. Horodecki, and 
K. Horodecki, Rev. Mod. Phys. 81, 865 (2009). 

[11] H. Li and F. D. M. Haldane, Phys. Rev. Lett. 101, 010504 
(2008). 


[12] A. J. A. James and R. M. Konik, Phys. Rev. B 87, 241103 
(2013). 

[13] V. Alba, M. Haque, and A. M. Lauchli, Phys. Rev. Lett. 
110, 260403 (2013). 

[14] M. B. Plenio and S. Virmani, Quant. Inf. Comput. 7, 1 
(2007). 

[15] O. F. Syljuasen, Phys. Lett. A 322, 25 (2004). 

[16] S.-J. Gu, G.-S. Tian, and H.-Q. Lin, Phys. Rev. A 71, 
052322 (2005). 

[17] S.-J. Gu, G.-S. Tian, and H.-Q. Lin, New J. Phys. 8, 61 
(2006). 

[18] C.-Y. Huang and F.-L. Lin, Phys. Rev. A 81, 032304 

( 2010 ). 

[19] J. Stasifiska, B. Rogers, M. Paternostro, G. De Ghiara, 
and A. Sanpera, Phys. Rev. A 89, 032330 (2014). 

[20] G. De Chiara, L. Lepori, M. Lewenstein, and A. Sanpera, 
Phys. Rev. Lett. 109, 237208 (2012). 

[21] J. 1. Girac and F. Verstraete, J. Phys. A: Math. Theor. 
42, 504004 (2009). 

[22] F. Verstraete, V. Murg, and J. Cirac, Adv. Phys. 57, 
143 (2008). 

[23] R. Orus, Ann. Phys. 349, 117 (2014). 

[24] S. R. White, Phys. Rev. Lett. 69, 2863 (1992). 





19 


[25] G. Vidal, Pliys. Rev. Lett. 91, 147902 (2003). 

[26] F. Verstraete and J. Cirac, (2004), 

arXiv:cond-mat/0407066. 

[27] A. Garcia-Saez and J. I. Latorre, Phys. Rev. B 87, 085130 
(2013). 

[28] H. G. Jiang, Z. Y. Weng, and T. Xiang, Phys. Rev. Lett. 
101, 090603 (2008). 

[29] H. H. Zhao, Z. Y. Xie, Q. N. Ghen, Z. G. Wei, J. W. Cai, 
and T. Xiang, Phys. Rev. B 81, 174411 (2010). 

[30] U. Schollwock, Ann. Phys. 326, 96 (2011). 

[31] Z.-G. Gu, M. Levin, and X.-G. Wen, Phys. Rev. B 78, 
205116 (2008). 

[32] J. Jordan, R. Oriis, G. Vidal, F. Verstraete, and J. I. 
Girac, Phys. Rev. Lett. 101, 250602 (2008). 

[33] R. Orus and G. Vidal, Phys. Rev. B 80, 094403 (2009). 

[34] S. Ostlund and S. Rommer, Phys. Rev. Lett. 75, 3537 
(1995). 

[35] S. Rommer and S. Ostlund, Phys. Rev. B 55, 2164 
(1997). 

[36] G. Vidal, Phys. Rev. Lett. 98, 070201 (2007). 

[37] M. Levin and G. P. Nave, Phys. Rev. Lett. 99, 120601 
(2007). 

[38] W. Magnus, Commun. Pure App. Math. 7, 649 (1954). 

[39] H. F. Trotter, Proc. Amer. Math. Soc. 10, 545 (1959). 

[40] R. Bamler, Diploma Thesis in Physics (2011), unpub¬ 
lished. 

[41] W. Li, J. von Delft, and T. Xiang, Phys. Rev. B 86, 
195137 (2012). 

[42] R. J. Baxter, J. Math. Phys. 9, 650 (1968). 

[43] T. Nishino and K. Okunishi, J. Phys. Soc. Jpn. 65, 891 
(1996). 

[44] T. Nishino and K. Okunishi, J. Phys. Soc. Jpn. 66, 3040 
(1997). 

[45] R. Orus, Phys. Rev. B 85, 205117 (2012). 

[46] R. Orris, Private communication. 

[47] B. Pirvu, V. Murg, J. 1. Girac, and F. Verstraete, New 
J. Phys. 12, 025012 (2010). 

[48] X.-N. Zhu and S.-M. Fei, Phys. Rev. A 90, 024304 (2014). 

[49] G. Gour, S. Bandyopadhyay, and B. C. Sanders, J. Math. 
Phys. 48, 012108 (2007). 

[50] P. Rungta, V. Buzek, C. M. Caves, M. Hillery, and G. J. 
Milburn, Phys. Rev. A 64, 042315 (2001). 

[51] B. Regula, S. Di Martino, S. Lee, and G. Adesso, Phys. 
Rev. Lett. 113, 110501 (2014). 

[52] E. Lieb, T. Schultz, and D. Mattis, Ann. Phys. 16, 407 
(1961). 

[53] H. W. J. Blote and Y. Deng, Phys. Rev. E 66, 066110 

( 2002 ). 

[54] F. Verstraete, K. Audenaert, J. Dehaene, and 
B. De Moor, J. Phys. A: Math. Gen. 34, 10327 (2001). 

[55] F. Verstraete, M. Popp, and J. 1. Cirac, Phys. Rev. Lett. 
92, 027901 (2004). 

[56] O. F. Syljuasen, Phys. Rev. A 68, 060301 (2003). 


[57] H.-J. Mikeska and A. K. Kolezhuk, in Quantum Mag¬ 
netism, Lecture Notes in Physics, Vol. 645, edited by 
U. Schollwock, J. Richter, D. J. Farnell, and R. F. Bishop 
(Springer Berlin Heidelberg, 2004) pp. 1-83. 

[58] L. Justino and T. R. de Oliveira, Phys. Rev. A 85, 052128 

( 2012 ). 

[59] V. S. Viswanath, S. Zhang, J. Stolze, and G. Muller, 
Phys. Rev. B 49, 9702 (1994). 

[60] S. Yunoki, Phys. Rev. B 65, 092402 (2002). 

[61] H.-Q. Lin, J. S. Flynn, and D. D. Betts, Phys. Rev. B 
64, 214411 (2001). 

[62] W.-L. You and Y.-L. Dong, Phys. Rev. B 84, 174426 

( 2011 ). 

[63] R. F. Bishop, D. J. J. Farnell, and J. B. Parkinson, J. 
Phys.: Cond. Matter 8, 11153 (1996). 

[64] C. N. Yang and C. P. Yang, Phys. Rev. 150, 321 (1966). 

[65] C. N. Yang and C. P. Yang, Phys. Rev. 150, 327 (1966). 

[66] C. N. Yang and C. P. Yang, Phys. Rev. 151, 258 (1966). 

[67] D. C. Mattis, The Theory of Magnetism I: Statics 

and Dynamics, Springer Series in Solid-State Sciences 
(Springer-Verlag, Berlin Heidelberg, 1981). 

[68] Z.-C. Gu and X.-G. Wen, Phys. Rev. B 80, 155131 
(2009). 

[69] A. W. Sandvik, Phys. Rev. B 56, 11678 (1997). 

[70] L.-A. Wu, M. S. Sarandy, and D. A. Lidar, Phys. Rev. 
Lett. 93, 250404 (2004). 

[71] Q. Meng and T. Dong-Ping, Chinese Phys. C 33, 249251 
(2009). 

[72] G.-H. Liu, W. Li, W.-L. You, G. Su, and G.-S. Tian, 
EPL 101, 57001 (2013). 

[73] S.-J. Gu, G.-S. Tian, and H.-Q. Lin, New J. Phys. 8, 61 
(2006). 

[74] D. Schwandt, F. Alet, and S. Capponi, Phys. Rev. Lett. 
103, 170501 (2009). 

[75] M. Jakob and J. A. Bergou, Opt. Comm. 283, 827 

( 2010 ). 

[76] M. A. Nielsen and 1. L. Chuang, Quantum Computation 
and Quantum Information (Cambridge University Press, 
Cambridge, England, 2000). 

[77] C. H. Bennett, D. P. DiVincenzo, J. A. Smolin, and 
W. K. Wootters, Phys. Rev. A 54, 3824 (1996). 

[78] S. Hill and W. K. Wootters, Phys. Rev. Lett. 78, 5022 
(1997). 

[79] W. K. Wootters, Phys. Rev. Lett. 80, 2245 (1998). 

[80] D. P. DiVincenzo, C. A. Fuchs, H. Mabuchi, J. A. Smolin, 
A. Thapliyal, and A. Uhlmann, LNCS 1509, 247 (1999). 

[81] T. Laustsen, F. Verstraete, and S. J. van Enk, Quant. 
Inf. Comp. 3, 64 (2003). 

[82] J. A. Smolin, F. Verstraete, and A. Winter, Phys. Rev. 
A 72, 052317 (2005). 

[83] G. Vidal and R. F. Werner, Phys. Rev. A 65, 032314 

( 2002 ). 

[84] M. B. Plenio, Phys. Rev. Lett. 95, 090503 (2005). 


